i = 0
gbm.results <- list()
ntree_seq <- seq(50,500,10)
for (year in start.year:end.year) {
  i = i + 1
  print(year)
  # Fit GBM
  gbm.results[[i]] <- list()
  dta.past <- dta[dta[,t.var]< year,]
  N <- nrow(dta.past)
  cv_results <- foreach (cv=1:cv.runs,
                         .combine='rbind',
                         .packages = c("gbm")) %dopar% {   # Loop over CV runs
                                         set.seed(52*cv)
                                         shuffle <- sample(1:N,
                                                           N,
                                                           replace=FALSE)
                                         predictions <- matrix(NA,nrow=N,
                                                               ncol=length(ntree_seq))
                                         sampleSplit = c(0,round((N/5*c(1:5))))
                                         for (f in 1:5) {
                                           # Fit lasso
                                           cv.fit = gbm.fit(x = as.matrix(dta.past[-shuffle[(sampleSplit[f]+1):(sampleSplit[f+1])],
                                                                                        rhs]),
                                                            y = as.matrix(dta.past[-shuffle[(sampleSplit[f]+1):(sampleSplit[f+1])],
                                                                                        dv]),
                                                            n.trees = max(ntree_seq),
                                                            distribution = gbm_count_params$distribution,
                                                            interaction.depth = gbm_count_params$interaction.depth,
                                                            n.minobsinnode = gbm_count_params$n.minobsinnode,
                                                            verbose = gbm_count_params$verbose,
                                                            shrinkage = gbm_count_params$shrinkage,
                                                            train.fraction = gbm_count_params$train.fraction,
                                                            bag.fraction = gbm_count_params$bag.fraction)
                                           predictions[shuffle[(sampleSplit[f]+1):(sampleSplit[f+1])],
                                                     ] <- predict(cv.fit,
                                                                  newdata = as.matrix(dta.past[shuffle[(sampleSplit[f]+1):(sampleSplit[f+1])],
                                                                                                    rhs]),
                                                                  type = "response",
                                                                  n.trees = ntree_seq)
                                         }
                                         cv_results <- c("CV Run"=0,
                                                         "Trees"=0,
                                                         "R2"=0)
                                         # Get error rates for each lambda used
                                         for (t in 1:length(ntree_seq)) {
                                           if (is.na(predictions[1,t])==0) {
                                             ###!####
                                             R2 <- Rdev(as.matrix(dta.past[,dv]),
                                                        predictions[,t])
                                             ### ! ####
                                             print(R2)
                                             if (R2>cv_results["R2"]) {
                                               cv_results <- c("CV Run"=cv,
                                                               "Trees" = ntree_seq[t],
                                                               "R2" = R2)
                                               best.t <- t
                                             }
                                           }
                                         }
                                         return(cv_results)
                                       }
  hours <- floor((proc.time()[3]-start.time)/3600)
  mins <- floor((proc.time()[3]-start.time)/60) - hours*60
  secs <- floor((proc.time()[3]-start.time)) - hours*3600 - mins*60
  
  print(paste(hours,
              "h",
              mins,
              "m",
              secs,
              "s elapsed",
              sep = ""))
  
  gbm.results[[i]]$trees <- round(mean(cv_results[,"Trees"]))
  
  mod = gbm.fit(x = as.matrix(dta[dta[,t.var]< year,
                                       rhs]),
                y = as.matrix(dta[dta[,t.var]< year,
                                       dv]),
                n.trees = gbm.results[[i]]$trees,
                distribution = gbm_count_params$distribution,
                interaction.depth = gbm_count_params$interaction.depth,
                n.minobsinnode = gbm_count_params$n.minobsinnode,
                verbose = gbm_count_params$verbose,
                shrinkage = gbm_count_params$shrinkage,
                train.fraction = gbm_count_params$train.fraction,
                bag.fraction = gbm_count_params$bag.fraction
  )
  gbm.results[[i]]$fit.oos <- predict(mod,
                                      newdata = as.matrix(dta[dta[,t.var] == year,
                                                                   rhs]),
                                      type = "response",
                                      n.trees = gbm.results[[i]]$trees)
  gbm.results[[i]]$actual.oos <- as.matrix(dta[dta[,t.var] == year,
                                                    dv])
  
  print(year)
}
save(gbm.results,
     file=paste(modeldir,"/",
                country,
                "_gbm_",
                "count",
                "_",
                rhs.group,fileext,
                ".RData",
                sep = ""))